#include <string>
#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "astro/cosmology.h"
#include "astro/cdmModel.h"
#include "astro/utilities.h"
#include <gsl/gsl_integration.h> 
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_rng.h>
#include <nr.h>
#include <nrD.h>
#include "time.h"

int const n=512;

/*
 * Cubic spline interpolation
 */
double InterpolateYvec(std:: vector<double> x, std:: vector<double> y,double xi){
  int n=x.size();
  if(x[n-1]>x[0]){
    if(xi>x[n-1]) return y[n-1];
    if(xi<x[0]) return y[0];
  }
  if(x[n-1]<x[0]){
    if(xi<x[n-1]) return y[n-1];
    if(xi>x[0]) return y[0];
  }
  int i = astro:: locate (x,xi);
  i = std::min (std::max (i,0), int (n)-2);
  double f=(xi-x[i])/(x[i+1]-x[i]);
  if(i>1 && i<n-2){
    double a0,a1,a2,a3,f2;
    f2 = f*f;
    a0 = y[i+2] - y[i+1] - y[i-1] + y[i];
    a1 = y[i-1] - y[i] - a0;
    a2 = y[i+1] - y[i-1];
    a3 = y[i];
    return a0*f*f2+a1*f2+a2*f+a3;
  }
  else return f*y[i+1]+(1-f)*y[i];
}

int main(){
  std:: string filein,str,fileout;
  std:: ifstream fin;
  filein = "INPUT";
  fin.open(filein.c_str());
  if(!fin){
    std:: cout << "INPUT file does not exist" << std:: endl;
    std:: cout << " I will STOP here!!!" << std:: endl;
    exit(1);
  }
  double wmin;
  double om0,oml,h0,sig8,z1,z2,w,lminmass,lmaxmass,howmanysky;
  int mftype,seedi;
  fin >> str;
  fin >> om0;
  fin >> str;
  fin >> oml;
  fin >> str;
  fin >> h0;
  fin >> str;
  fin >> sig8;
  fin >> str;
  fin >> w;
  fin >> str;
  fin >> z1;
  fin >> str;
  fin >> z2;
  fin >> str;
  fin >> lminmass;
  fin >> str;
  fin >> lmaxmass;
  fin >> str;
  fin >> fileout;
  fin >> str;
  fin >> howmanysky;
  fin >> str;
  fin >> mftype;
  fin >> str;
  fin >> wmin;
  fin >> str;
  fin >> seedi;
  fin.close();  

  if(seedi==0){
    seedi = time(NULL);
    std:: cout << " random seed associated to the time " << seedi << std:: endl;
  }

  if(wmin<lminmass || wmin>lmaxmass){
    std:: cout << " mass to write in the output file is not between " << std:: endl;
    std:: cout << " minMass and maxMass " << std:: endl;
    std:: cout << " mass to write = " << pow(10.,wmin) << std:: endl;
    std:: cout << " minMass = " << pow(10.,lminmass) << std:: endl;
    std:: cout << " maxMass = " << pow(10.,lmaxmass) << std:: endl;
    std:: cout << " check it out, I will STOP here !!!" << std:: endl;
    exit(1);
  }

  // read in put file end initialise the cosmological model and vectors
  astro:: cosmology co(om0,oml,h0);
  astro:: cdmModel cdm(co,8.,sig8);  
  std:: cout << " cosmology: " << std:: endl;
  std:: cout << "     om0="<<om0<<"   oml="<<oml<< std::endl;
  std:: cout << "     h0="<<h0<<"   sig8="<<sig8<< "   w="<<w<< std::endl;
  std:: cout << "   " << std:: endl;
  std:: cout << "   generating list of haloes from the " << std:: endl;
  if(mftype == 1) std:: cout << " Sheth and Tormen mass function " << std:: endl;
  if(mftype == 0) std:: cout << " Press and Shechter mass function " << std:: endl;
  if(mftype == 2) std:: cout << " Jenkins et al. mass function " << std:: endl;
  std:: cout << " in the redshift bin = " << z1 << " - " << z2 << std:: endl;
  std:: cout << " writing them in " << fileout << std:: endl;
 
  double tiny = 1.e-4;
  std :: vector<double> logm(n),m(n),distM(n),logdistM(n);
  astro:: fill_linear(logm,n,lminmass,lmaxmass);
  for(int i=0;i<n;i++ ) m[i] = pow(10.,logm[i]);
  double mMin = pow(10.,lminmass);
  double normalization = cdm.number(mMin,z1,z2,mftype);
  std ::  ofstream Diff;
  Diff.open("distM.dat");
  for (int i=0; i<n; i++){
    distM[i] = cdm.number(m[i],z1,z2,mftype)/normalization;
    Diff << distM[i] << "   " << logm[i] << "  "   << cdm.number(m[i],z1,z2,mftype) 
	 << std :: endl;
  }
  std:: vector<double> logm2,distM2;
  for(int i=0;i<n-1;i++){
    if(i==0){
      distM2.push_back((distM[i]));
      logm2.push_back(logm[i]);
    }
    else{
      if((distM[i]-distM[i-1])>tiny){
	distM2.push_back((distM[i]));
	logm2.push_back(logm[i]);
      }
    }
  }

  for(int i=0;i<n;i++){
    if(distM[i]<=0){
      std:: cout<< "distM  = "<< distM[i] << "for logM " << logm[i] << std:: endl;
      exit(1);
    }
    logdistM[i] = log10(distM[i]);
  }

  long seed = long(seedi);
  //random number and random number of haloes
  const gsl_rng_type *T;
  gsl_rng *r;
  double Nd;
  std:: cout << lminmass << std:: endl;
  Nd = cdm.number(pow(10.,lminmass),z1,z2,mftype);
  std:: cout << Nd << std:: endl;
  Nd = Nd*howmanysky;
  long N = long (poidev(long(Nd), &seed) + 0.5);
  std:: cout << N << std:: endl;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);  
  gsl_rng_set(r,int(2*seedi*seedi));           
  int i = 0;
  std ::  ofstream Masses;
  Masses.open(fileout.c_str());
  int ji = 0;
  while (i<N){
    double u = gsl_rng_uniform (r);
    if(u>0.){
      u = log10(u);
      double v = InterpolateYvec(logdistM,logm,u);
      if(v>=wmin){
	Masses << pow(10.,v) << std :: endl;
	ji++;
      }
      i++;
    }
  }    
  Masses.close();
  std:: cout << ji << std:: endl;
  gsl_rng_free (r);
}
